 function [GG,RR,C,eu,SDX,ZZ,initss,ssvec,flag,ssnames,stanames,shonames]=modchimodDtFacpInfExpGapS9(parest,~,~,calval,posest,poscal)
% =========================================================================
%
% Taken from modchimodDtFacpInfExpGap_obs
% This code is trimmed to 9 signals 
% and requires to pass the parameters of the signals as specific
% coefficients 
%
% ADDED additional parameter that determines how many policy signals are
% observed 
%
% Checked identical to code used for estimation Sep 7 2011 
% =========================================================================

if nargout < 11 
    flag_repstanames=0; 
else 
    flag_repstanames=1; 
end 
if nargout < 10 
    flag_repssnames=0; 
else
    flag_repssnames=1; 
end

flag.ssok=1;
flag.solver=1; 

initss=[]; 
ssnames=[]; 
ssvec  =[]; 


% if nargin~=1 && nargin~=4 
%     error('Need either 1 or 4 inputs') 
% end 

% =========================================================================
%% *PART 1. Declare All Variables in the Model Solution $$ S_{t} $$* 

%__________________________________________________________________________
% 1.a Endogenous Variables 
y=1;            % output
k=2;            % capital
L=3;            % hours
Rk=4;           % return on capital
w=5;            % real wages
x=6;            % marginal utility of labor
p=7;            % inflation
s=8;            % marginal cost
lambda=9;       % multiplier
c=10;           % consumption
R=11;           % interst rate
u=12;           % capital utilization
phi=13;         % multiplier
i=14;           % investment
kbar=15;        % "gross" capital
%__________________________________________________________________________
% 1.b Exogenous driving processes
mp=16;
z =17;          % Productivity shock
g =18;          % Public spending shock
miu=19;         % Relative price (non stationary), this is Growth Rate
lambdap=20;     % Markup shock AR(1) 
psi=21;         % Preference shock AR(1) 
b=22;           % Discount factor shock 
upsilon=23;     % MEI, stationary, investment shock
pitarg=24;      % Inflation Drift
wmarkwn=25;     % White Noise Wage Markup 
%__________________________________________________________________________
% 1.c Observable growth rates (nominal) 
ynomobs=26; 
cnomobs=27;
inomobs=28; 
wnomobs=29; 
%__________________________________________________________________________
% 1.c.a Auxiliary Variables or Observable Prices
yrobs   =30; 
gdp     =31; 
dpinv   =32; 
gdpdef  =33; 
crobs   =34;   % delta(Nominal C - Pc)
irobs   =35;   % delta(Nominal I - Pc)
wrobs   =36;   % delta(Nominal W - Pc)
ygap    =37;   % Output Gap 
exreal1 =38;   % Expected Real Rate next period 
%__________________________________________________________________________
% 1.d Lags
yrobs_1=39;   % Real Observable Output (adjusted for utilization) lag1      
yrobs_2=40;   % Real Observable Output (adjusted for utilization) lag2
yobsan =41;   % Annualize Observable Output 
p_1    =42;   % Lagged C deflator 
p_2    =43;   % Lagged C deflator 
dpma   =44;   % Four quarter inflation  
dagtrend=45;  % First difference of the aggregate trend 
% Make DAGTREND Last Variable here for automatic updating of the next
% position 
%__________________________________________________________________________
% 1.e 1 through 12 step ahead expectations of FFR 
startRlead=dagtrend+1;  % Number to start with CURRENTLY 46
R_lead1=startRlead;  
R_lead2=startRlead+1; 
R_lead3=startRlead+2; 
R_lead4=startRlead+3;
R_lead5=startRlead+4;
R_lead6=startRlead+5;
R_lead7=startRlead+6;
R_lead8=startRlead+7;
R_lead9=startRlead+8;
R_lead10=startRlead+9;
R_lead11=startRlead+10;
R_lead12=startRlead+11;
Rlead_vec=(R_lead1:R_lead12); 
%__________________________________________________________________________
% 1.f Policy signals and Inflation Expectations
startSig=startRlead+12; % Number to start with CURRENTLY 58 
sigR_0=startSig; 
sigR_1=startSig+1; 
sigR_2=startSig+2; 
sigR_3=startSig+3; 
sigR_4=startSig+4; 
sigR_5=startSig+5; 
sigR_6=startSig+6; 
sigR_7=startSig+7; 
sigR_8=startSig+8; 
sigR_9=startSig+9; 

posExpInf = startSig+9+1;
pe_vec = (posExpInf):(posExpInf+39);
expectedPi40 = posExpInf+40;
realRate12q  = posExpInf+41; 
%__________________________________________________________________________
% 1.g Flex Price and Wage Economy * 
startstar=posExpInf+41+1;        % Number to start with CURRENTLY 111
ystar=startstar;              % output
kstar=startstar+1;            % capital
Lstar=startstar+2;            % hours
Rkstar=startstar+3;           % return on capital
wstar=startstar+4;            % real wages
xstar=startstar+5;            % marginal utility of labor
Rstar=startstar+6;            % interst rate
sstar=startstar+7;            % marginal cost
lambdastar=startstar+8;       % multiplier
cstar=startstar+9;            % consumption
ustar=startstar+10;           % capital utilization
phistar=startstar+11;         % multiplier
istar=startstar+12;           % investment
kbarstar=startstar+13;        % "gross" capital
gdpstar=startstar+14;         % GDP 
NY=startstar+14;

% =========================================================================

% =========================================================================
%% *PART 2. Declare Innovations to the exogenous variables $$ \eta_{t} $$* 

nsignals = 9;    % Number of policy signals 
NX=10+ nsignals; 
%__________________________________________________________________________
% 2.a Structural Disturbances 
Rs=1;           % MP
zs=2;           % Technology 
gs=3;           % Public spending
mius=4;         % Investment     
lambdaps=5;     % Mark-up
psis=6;         % Leisure preference
bs=7;           % Discount factor
upsilons=8;     % investment shock
shpitarg=9;     % Inflation Drify
wmarkwnsh=10; 
%__________________________________________________________________________
% 2.b Policy Signals 
Rs_lag1=11;    % 1 period ahead
Rs_lag2=12;    % 2 periods ahead
Rs_lag3=13;    % 3 periods ahead
Rs_lag4=14;    % 4 periods ahead
Rs_lag5=15;    % 5 periods ahead
Rs_lag6=16;    % 6 periods ahead
Rs_lag7=17;    % 7 periods ahead
Rs_lag8=18;    % 8 periods ahead
Rs_lag9=19;    % 9 periods ahead
% =========================================================================

% =========================================================================
%% *PART 3. Declare Parameters* 
ncof=        68;    % Number of coefficients 
numpar=      ncof;  % Equal to Parameters, distinction irrelevant here 
if nargin < 4 || isempty( calval ) 
    param=parest  ; 
else
    param=zeros(numpar,1); 
    param(posest)=parest; 
    param(poscal)=calval; 
end; 
%__________________________________________________________________________
% 3.1 Parameters for all but multiple prices, GDP and I deflators 

alpha=param(1);         % share of capital in the prod. function
delta=param(2);         % capital depreciation rate
iotap=param(3);         % price indexation (=0 is static indexation, =1 is dynamic)
iotaw=param(4);         % wages indexation
gammastar100=param(5);  % steady state growth rate of technology
gammamiu100=param(6);   % steady state growth rate of capital embodied technology
h=param(7);             % habit formation parameter
lambdapss=param(8);     % steady state of mark-up shock (pins down steady state of wages)
lambdaw=param(9);       % wage mark-up
Lss=param(10);          % steady state for leisure (the ss for leisure is pinned down by psiss. convenient to parameterize the model in terms of Lss)
pss100=param(11);       % steady state quarterly inflation (multiplied by 100)
Fbeta=param(12);        % weird parameter of SW that is a function of beta and the steady state quarterly real rate of interest (ensures that beta is less than 1)
gss=param(13);          % steady state of the public spending shock
niu=param(14);          % curvature of the utility function for leisure
xip=param(15);          % price stickiness
xiw=param(16);          % wage stickiness
chi=param(17);          % elasticity of the capital utilization cost function
S=param(18);            % investment adjustment cost
fp=param(19);           % reaction to annualized inflation
fy=param(20);           % reaction to annualized output growth 
rhoR=param(21);         % policy smoothing  
% AR coefficients of exogenous prorcesses 
rhoz=param(22); 
rhog=param(23); 
rhomiu=param(24); 
rholambdap=param(25); 
rhopsi=param(26); 
rhob=param(27); 
rhopitarg=param(28); 
rhomp=param(29); 
rhoupsilon=param(30); 
% Std of structural shocks 
sdstruct=param(31:40);
kappa       = -param(41);   % Scaling for spredd  
sdspreadme  =  param(42);   % Measurement error spread
% Std Policy Signals
% ===============================================================
% KEYWORD: NewWithS9
%% Requires std to the 9 policy signals to be defined as parameters 
stdsignals=param(43:51); 
if length(stdsignals)~=nsignals; error('Wrong Number of Signals');end 

% Number of ExFFRObserved 
NExFFRObserved=param(52); 

%% Parameters for inflation expectations 
posCPI=param(53); 
expInfLoading=param(54); 
expInfCons   =param(55); 
end_nonp=55; 
% _________________________________________________________________________
% 3.2  Begin Parameters for Multiple Prices 
%__________________________________________________________________________
% 3.2.1  How to aggregate gdp deflator
flag_wadd1  =param(end_nonp+1); 
nserp       =param(end_nonp+2); 
% __________________________________________________________________
% 3.2.1 Declare positions for 
% Rows for 
% a) loadings 
% b) ar ID errors 
% c) constants 
% d) volatilities 
%
% *Note:* There are NSERP+1 loadings. First 2 are for the GDP Deflator.
% 1st of W(1)*Price consumption.2nd on W(2)*Prince Investment. 
% Remaining are for other prices on the Consumption Deflators 
row_load=(end_nonp+3):(end_nonp+nserp+3); 
row_ar  =row_load(end)*ones(1,nserp)+(1:nserp); 
row_std =row_ar(end)*ones(1,nserp)+(1:nserp);
row_c   =row_std(end)*ones(1,nserp)+(1:nserp); 

%__________________________________________________________________________
% 3.2.3 vector of factor loadings 
bet_gdpdef  =param(row_load(1:2)); 
bet_c       =param(row_load(3:end)); 

bet_CPI     =bet_c(posCPI-1); 
if abs(expInfLoading-999) < 0.01; 
    expInfLoading=bet_CPI;
end



%__________________________________________________________________________
% 3.2.4 vector of idiosyncratic errors 
rhoid=2*param(row_ar)-1; 
if any(rhoid < -0.995)~=0 || any(rhoid > 0.995)~= 0 
    disp('Explosive AR ID roots'); 
    GG=[];C=[];RR=[];eu=[0;0];SDX=[];ZZ=[];
    return
end

%__________________________________________________________________________
% 3.2.5 vector of constants 
cprices =param(row_c); 
if all( cprices ~= 0 ) 
    error('At least one of the constants must be set equal to 0')
end
constantCPI=cprices(posCPI); 

%__________________________________________________________________________
% 3.2.6 vector of idiosyncratic errors
sdid             =param(row_std); 

%__________________________________________________________________________
% 3.2.7 constants for I deflator 
cons_idef        =param(row_c(end)+1:end);
if length(cons_idef)> 1 
    error('Wrong Dimension of Investment Deflator') 
end 

%==========================================================================
%% *PART 4. Compute Steady State (SS)* 
 
gammastar=gammastar100/100;
gammamiu=gammamiu100/100;
beta=100/(Fbeta+100);
rss=exp(gammastar)/beta-1;
rss100=rss*100;
gss=1/(1-gss);
Rkss=exp(gammastar+gammamiu)/beta-1+delta;
sss=1/(1+lambdapss);
wss=(sss*((1-alpha)^(1-alpha))/((alpha^(-alpha))*Rkss^alpha))^(1/(1-alpha));
kLss=(wss/Rkss)*alpha/(1-alpha);
FLss=(kLss^alpha-Rkss*kLss-wss);
yLss=kLss^alpha-FLss;

% Scale of the economy is fixed at 1, i.e. independent of LSS which now
% only enters observation equation for hours 
Lsscale=1;
kss=kLss*Lsscale;
iss=(1-(1-delta)*exp(-gammastar-gammamiu))*kss*exp(gammastar+gammamiu);
F=FLss*Lsscale;
yss=yLss*Lsscale;
css=yss/gss-iss;

ssvec  =[css/yss;iss/yss;kss/yss;rss100;100*Rkss]; 
if  flag_repssnames==1;
    ssnames={'c/y'  ;'i/y'  ;'k/y';  'rss100';'Rkss100'};
end
    
%%==========================================================================

%==========================================================================
%% *PART 5. Declare System Matrices* 
% Model to be solved is written as 
% $$ A_{0}Y_{t+1} + A_{1}Y_{t}=A_{2}Y_{t-1} + PSI \eta_{t} $$  
A0  = zeros(NY,NY) ;
A1  = zeros(NY,NY) ;
A2  = zeros(NY,NY) ;
PSI = zeros(NY,NX) ;
%==========================================================================
expg=exp(gammastar);

%% *PART 6. Begin Filling Equation Rows* 

%% Eq 1, Output (y) from production function 
A1(y,y)=1;
A1(y,k)=-((yss+F)/yss)*alpha;
A1(y,L)=-((yss+F)/yss)*(1-alpha);

%% Eq 1*, Output (y*) in flex price & wage 
A1(ystar,ystar)=1;
A1(ystar,kstar)=-((yss+F)/yss)*alpha;
A1(ystar,Lstar)=-((yss+F)/yss)*(1-alpha);

%% Eq. 2, Effective Capital (K) 
A1(k,k)=1;
A1(k,u)=-1;
A1(k,z)=1;
A2(k,kbar)=1;
A1(k,miu)=alpha/(1-alpha)+1;

%% Eq 2*, Effective Capital (K*) in flex price & wage 
A1(kstar,kstar)=1;
A1(kstar,ustar)=-1;
A1(kstar,z)=1;
A2(kstar,kbarstar)=1;
A1(kstar,miu)=alpha/(1-alpha)+1;

%% Eq. 3 Hours (L), from cost minimization 
A1(L,Rk)=1;
A1(L,w)=-1;
A1(L,k)=1;
A1(L,L)=-1;

%% Eq. 3* Hours (L*) in flex price & wage 
A1(Lstar,Rkstar)=1;
A1(Lstar,wstar)=-1;
A1(Lstar,kstar)=1;
A1(Lstar,Lstar)=-1;

%% Eq. 4 Return to Capital (Rk) from utilization FOC
A1(Rk,Rk)=1;
A1(Rk,u)=-chi;

%% Eq. 4* Return to Capital (Rk*) in flex price & wage
A1(Rkstar,Rkstar)=1;
A1(Rkstar,ustar)=-chi;

%% Eq. 5, Real Wages (w) PC (with split of shocks)
slopewages=(1-beta*xiw)*(1-xiw)/((1+beta)*xiw*(1+niu*(1+1/lambdaw)));
A1(w,w)=1;
A0(w,w)= -beta/(1+beta);
A0(w,p)= -beta/(1+beta);
A1(w,x)=slopewages;
A1(w,z)=(1+iotaw*beta-rhoz*beta)/(1+beta);
A1(w,miu) =(alpha/(1-alpha))*((1+iotaw*beta-rhomiu*beta)) *(1/(1+beta)); 

A1(w,p)=(1+iotaw*beta)/(1+beta);
A2(w,p)=iotaw/(1+beta);
A2(w,w)=1/(1+beta);

A2(w,z)=iotaw/(1+beta);
A2(w,miu)=(alpha*iotaw)/( (1+beta)*(1-alpha) ); 

A1(w,wmarkwn)= -1; 

%% Eq. 5*, Real Wages (w*) in flex price & wage 
% $$ x^{*}_{t} = 0 $$ i.e. wage gap is zero, so workers are always on their supply curve 
A1(wstar,xstar)=1;

%% Eq. 6, Wage Gap (x), (MRS-Real Wage)
A1(x,x)=1;
A1(x,w)=-1;
A1(x,lambda)=-1;
A1(x,b)=1/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2));
A1(x,L)=niu;
A1(x,psi)=1; 

%% Eq. 6*, Wage Gap (x*), (MRS-Real Wage) in flex price & wage 
% Equal zero from eq 5*
A1(xstar,xstar)=1;
A1(xstar,wstar)=-1;
A1(xstar,lambdastar)=-1;
A1(xstar,b)=1/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2));
A1(xstar,Lstar)=niu;
A1(xstar,psi)=1; 

%% Eq 7, Prices (p) PC 
A1(p,p)=1;
A0(p,p)=-beta/(1+iotap*beta);
A2(p,p)=iotap/(1+iotap*beta);
A1(p,s)=-((1-beta*xip)*(1-xip)/((1+iotap*beta)*xip));
A1(p,lambdap)=-1;       

%% Eq 7*, Prices (p*) in flex price & wage
% $$ s^{*}_{t} = 0 $$ i.e. deviations of marginal cost from SS are zero, so
% firms are always on their supply curve
A1(Rstar,sstar)=1;

%% Eq 8, Marginal Cost (s)
A1(s,s)=1;
A1(s,Rk)=-alpha;
A1(s,w)=-(1-alpha);

%% Eq 8*, Marginal Cost (s*) in flex price & wage 
A1(sstar,sstar)=1;
A1(sstar,Rkstar)=-alpha;
A1(sstar,wstar)=-(1-alpha);

%% Eq 9, Lagrange Multiplier IBC (Lambda) 
A1(lambda,lambda)=1;
A1(lambda,R)=-1;
A0(lambda,lambda)=-1;
A0(lambda,p)=1;
A1(lambda,z)=rhoz;
A1(lambda,miu)=(rhomiu)*alpha/(1-alpha);

%% Eq 9*, Lagrange Multiplier (Lambda*) in flex price & wage
A1(lambdastar,lambdastar)=1;
A1(lambdastar,Rstar)=-1;
A0(lambdastar,lambdastar)=-1;
A1(lambdastar,z)=rhoz;
A1(lambdastar,miu)=(rhomiu)*alpha/(1-alpha);

%% Eq 10, Consumption (c)
A1(c,lambda)=(expg-h*beta)*(expg-h);
A1(c,b)=-(expg-h*beta*rhob)*(expg-h)/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2)); 
A1(c,z)=-(beta*h*expg*rhoz-h*expg);
A1(c,miu)=(-(beta*h*expg*rhomiu-h*expg))*alpha/(1-alpha);
A1(c,c)=(expg^2+beta*h^2);
A0(c,c)=-beta*h*expg;
A2(c,c)=expg*h;

%% Eq 10*, Consumption (c*) in flex price & wage
A1(cstar,lambdastar)=(expg-h*beta)*(expg-h);
A1(cstar,b)=-(expg-h*beta*rhob)*(expg-h)/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2));  
A1(cstar,z)=-(beta*h*expg*rhoz-h*expg);
A1(cstar,miu)=(-(beta*h*expg*rhomiu-h*expg))*alpha/(1-alpha);
A1(cstar,cstar)=(expg^2+beta*h^2);
A0(cstar,cstar)=-beta*h*expg;
A2(cstar,cstar)=expg*h;

%% Eq 11, Taylor Rule for Nominal Interest Rate (R)
A1(R,R)=1;
% 10.1 Systematic Component 
A2(R,R)=rhoR;
A1(R,dpma)  = -(1-rhoR)*fp; % Annualized Inflation 
A1(R,gdp )  = -(1-rhoR)*fy; % Annualized Observable Output 
% 10.2 Non-Systematic Component 
A1(R,mp)      = -1;         % Contemporaneous Innovation 
A1(R,sigR_0)  = -1;         % Signals 
A1(R,pitarg)  =(1-rhoR)*fp; % Inflation Drift

%% Eq. 12, Utilization (U) from market clearing condition 
A1(u,c)=css/yss;
A1(u,i)=iss/yss;
A1(u,y)=-1/gss;
A1(u,g)=1/gss;
A1(u,u)=kss*Rkss/yss;

%% Eq. 12*, Utilization (U) in flex price and wage 
A1(ustar,cstar)=css/yss;
A1(ustar,istar)=iss/yss;
A1(ustar,ystar)=-1/gss;
A1(ustar,g)=1/gss;
A1(ustar,ustar)=kss*Rkss/yss;

%% Eq. 13, Shadow Value of Capital (phi)
A1(phi,phi)=1;
A0(phi,phi)=-beta*exp(-gammastar-gammamiu)*(1-delta);
A1(phi,z)=rhoz;
A1(phi,miu)=(rhomiu)*(alpha/(1-alpha)+1);
A0(phi,lambda)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));
A0(phi,Rk)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));

%% Eq. 13*, Shadow Value of Capital (phi*) in flex price and wage
A1(phistar,phistar)=1;
A0(phistar,phistar)=-beta*exp(-gammastar-gammamiu)*(1-delta);
A1(phistar,z)=rhoz;
A1(phistar,miu)=(rhomiu)*(alpha/(1-alpha)+1);
A0(phistar,lambdastar)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));
A0(phistar,Rkstar)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));

%% Eq. 14, Investment (i, in consumption units) 
expgmiu=exp(gammastar+gammamiu);
A1(i,lambda)=1/(S*expgmiu^2);
A1(i,phi)=-1/(S*expgmiu^2);
A1(i,upsilon)=-1/(S*expgmiu^2);
A1(i,miu)=((1-beta*rhomiu))*(alpha/(1-alpha)+1); %no normalization
A1(i,i)=(1+beta);
A1(i,z)=(1-beta*rhoz);
A0(i,i)=-beta;
A2(i,i)=1;

%% Eq. 14*, Investment (i*, in consumption units) in flex price and wage 
A1(istar,lambdastar)=1/(S*expgmiu^2);
A1(istar,phistar)=-1/(S*expgmiu^2);
A1(istar,upsilon)=-1/(S*expgmiu^2);
A1(istar,miu)=((1-beta*rhomiu))*(alpha/(1-alpha)+1); %no normalization
A1(istar,istar)=(1+beta);
A1(istar,z)=(1-beta*rhoz);
A0(istar,istar)=-beta;
A2(istar,istar)=1;

%% Eq 15, Capital accumulation (kbar)
% Note: Recall, effective capital above is KBAR*UTILIZATION 
A1(kbar,kbar)=1;
A1(kbar,i)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A1(kbar,upsilon)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A2(kbar,kbar)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbar,z)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbar,miu)=((1-delta)*exp(-gammastar-gammamiu))*(alpha/(1-alpha)+1);

%% Eq 15*, Capital accumulation (kbar*) in flex price and wage 
A1(kbarstar,kbarstar)=1;
A1(kbarstar,istar)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A1(kbarstar,upsilon)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A2(kbarstar,kbarstar)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbarstar,z)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbarstar,miu)=((1-delta)*exp(-gammastar-gammamiu))*(alpha/(1-alpha)+1);

%% Eq 16, MP Exogenous Policy Innovation (unexpected) 
A1(mp,mp)=1; A2(mp,mp)=rhomp; PSI(mp,Rs)=1; 

%% Eq 17, Exogenous Growth Rate in Neutral Technoloy 
A1(z,z)=1; A2(z,z)=rhoz; PSI(z,zs)=1;

%% Eq 18, Exogenous Process in G Spending 
A1(g,g)=1; A2(g,g)=rhog; PSI(g,gs)=1;

%% Eq 19, Exogenous Process in Investment Specific Tech or Relative Price C
% to I 
A1(miu,miu)=1; A2(miu,miu)=rhomiu; PSI(miu,mius)=1;

%% Eq 20, Exogenous Process in Price Markup 
A1(lambdap,lambdap)=1; A2(lambdap,lambdap)=rholambdap;PSI(lambdap,lambdaps)=1;

%% Eq 21, Exogenous Process in Labor Disutility (AR1) 
A1(psi,psi)=1; A2(psi,psi)=rhopsi; PSI(psi,psis)=1;

%% Eq 22, Exogenous Process in Intertemporal Preference
A1(b,b)=1; A2(b,b)=rhob; PSI(b,bs)=1;

%% Eq 23, Exogenous Process in MEI 
A1(upsilon,upsilon)=1; A2(upsilon,upsilon)=rhoupsilon; PSI(upsilon,upsilons)=1;

%% Eq 24, Exogenous Process in Inflation Drift 
A1(pitarg,pitarg)=1; A2(pitarg,pitarg)=rhopitarg; PSI(pitarg,shpitarg)=1;

%% Eq 25, Exogenous Process in Wage Markup (white noise) 
A1(wmarkwn,wmarkwn)=1; PSI(wmarkwn,wmarkwnsh)=1; 

%% Eq. 26, Observable GDP Nominal (add P back)
A1(ynomobs,ynomobs) =  1;
A1(ynomobs,gdp)     = -1; 
A1(ynomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(ynomobs,p)       = -1; 

A2(ynomobs,gdp)     = -1; 

%% Eq. 27, Observable Consumption Nominal (add P back)
A1(cnomobs,cnomobs) =1;
A1(cnomobs,c)       = -1; 
A1(cnomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(cnomobs,p)= -1; 

A2(cnomobs,c)= -1; 

%% Eq. 28 Observable Investment Nominal (add P back)
% In consumption units
A1(inomobs,inomobs)   =1;
A1(inomobs,i)         = -1; 
A1(inomobs,dagtrend)  = -1; 
% Add inflation, since nominal 
A1(inomobs,p)= -1; 

A2(inomobs,i)= -1; 

%% Eq. 29 Observable Hourly Wages  Nominal (add P back)
A1(wnomobs,wnomobs) =1;
A1(wnomobs,w)       = -1; 
A1(wnomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(wnomobs,p)= -1; 

A2(wnomobs,w)= -1; 

%% Eq. 30 Observable Real GDP Growth (in C Units)
% Corrected observation equation for utilization 
% $$ yrobs_{t} = gdp_{t}-gdp_{t-1} + agtrend_{t} $$ 
% where AGTREND is the aggregate trend 
A1(yrobs,yrobs)=1; 
A1(yrobs,gdp)= -1;
A1(yrobs,dagtrend)= -1; 
A2(yrobs,gdp)= -1; 

%% Eq. 31 GDP (detrended), Real Detrended Output corrected for Utilization 
A1(gdp,gdp)= -1; 
A1(gdp,y)=1; 
A1(gdp,u)= -kss*Rkss/yss; 

%% Eq. 31* GDP* in flex price and wage
A1(gdpstar,gdpstar)= -1; 
A1(gdpstar,ystar)=1; 
A1(gdpstar,ustar)= -kss*Rkss/yss; 

%% Eq. 32 Growth Rate in Investment Deflator 
% $$ \Delta P^{I}_{t} = \Delta P^{C}_{t} - \mu_{t} $$
A1(dpinv,dpinv)= -1; 
A1(dpinv,p)    =  1; 
A1(dpinv,miu)  = -1; 

%% Eq. 33 GDP Deflator 
% $$ P^{GDP}_{t} = \beta_{1} P^{C}_{t} + \beta_{2} P^{I}_{t} $$ 
% 
% Recall there are 2 ways to define the weights 
%
% *First:* $$ {\frac{C}{C+I} , \frac{I}{C+I} } $$ 
% 
% *Second:* $$ {\frac{C}{Y}  , \frac{I}{Y}   } $$
%
% Note that in second case weights will not add up to one, which is not
% desirable, but opens up the door for adding a 3 deflator, i.e. price of G
% and NX. 
bet_gdpdef=bet_gdpdef(:); 
weights=zeros(1,2); 
if flag_wadd1==1
    weights(1)=css/(css+iss);
    weights(2)=iss/(css+iss);
else
    weights(1)=css/yss;
    weights(2)=iss/yss;
end
A1(gdpdef,gdpdef) = -1; 
A1(gdpdef,p)      = weights(1)*bet_gdpdef(1); 
A1(gdpdef,dpinv)  = weights(2)*bet_gdpdef(2); 


%% Eq. 34, Nominal C deflated by Model-based Pc deflator 
A1(crobs,crobs)  =-1; 
A1(crobs,cnomobs)=1; 
A1(crobs,p      )= -1; 

%% Eq. 35, Nominal I deflated by Model-based Pc deflator 
A1(irobs,irobs)  =-1; 
A1(irobs,inomobs)=1; 
A1(irobs,p      )= -1; 

%% Eq. 36, Nominal I deflated by Model-based Pc deflator 
A1(wrobs,wrobs)  =-1; 
A1(wrobs,wnomobs)=1; 
A1(wrobs,p      )= -1; 

%% Eq. 37, Output gap $$ GDP_{t} - GDP^{*}_{t} $$ 
A1(ygap,ygap)     =-1; 
A1(ygap,gdp)      =1; 
A1(ygap,gdpstar)  = -1; 
 
%% Eq 38 EXREAL1: Expected Real Rate next period $$ E[ R_{t} - \pi_{t+1} | t ] $$
A1(exreal1,exreal1) = -1; 
A1(exreal1,R)       =1; 
A0(exreal1,p)       = -1; 

%% Eq. 39 Lag 1 Observable Real GDP Growth (in C Units)
A1(yrobs_1,yrobs_1)=1; A2(yrobs_1,yrobs)=1;

%% Eq. 40 Lag 2 Observable Real GDP Growth (in C Units)
A1(yrobs_2,yrobs_2)=1; A2(yrobs_2,yrobs_1)=1;

%% Eq. 41 Annualized Observable Real GDP Growth (in C Units)
% $$ \Delta y^{A}_{t}= \sum_{j=0}^{3}yrobs_{t-j} $$
A1(yobsan,yobsan)= -1;

A1(yobsan,yrobs)  =1; 
A1(yobsan,yrobs_1)=1; 
A1(yobsan,yrobs_2)=1;
A2(yobsan,yrobs_2)= -1;

%% Eq. 42 Lag 1 Inflation (C)
A1(p_1,p_1)=1; A2(p_1,p)=1;
%% Eq. 43 Lag 2 Inflation (C)
A1(p_2,p_2)=1; A2(p_2,p_1)=1;
%% Eq. 44 Annualized Inflation (C) 
% $$ \Delta P^{A,C}_{t}= \frac{ \sum_{j=0}^{3} \Delta P^{C}_{t-j} }{4} $$
A1(dpma,dpma)= -1;
A1(dpma,p)   = 1/4;
A1(dpma,p_1) = 1/4;
A1(dpma,p_2) = 1/4;
A2(dpma,p_2) = -1/4;
%% Eq. 45, Growth Rate in Aggregate Trend 
% $$ \Gamma_{t} = Z_{t} + \frac{\alpha}{1-\alpha} \mu_{t} $$ 
A1(dagtrend,dagtrend)= -1; 
A1(dagtrend,z)=1; 
A1(dagtrend,miu)=alpha/(1-alpha); 

%% Eq. 46-58 Expectations of Future Interest Rates 
% $$ E[ R_{t+j} | t ] j=1,..,4 $$ In that order 
A1(R_lead1,R_lead1)=1; A0(R_lead1,R)= -1; 
A1(R_lead2,R_lead2)=1; A0(R_lead2,R_lead1)= -1;
A1(R_lead3,R_lead3)=1; A0(R_lead3,R_lead2)= -1;
A1(R_lead4,R_lead4)=1; A0(R_lead4,R_lead3)= -1;
A1(R_lead5,R_lead5)=1; A0(R_lead5,R_lead4)= -1;
A1(R_lead6,R_lead6)=1; A0(R_lead6,R_lead5)= -1;
A1(R_lead7,R_lead7)=1; A0(R_lead7,R_lead6)= -1;
A1(R_lead8,R_lead8)=1; A0(R_lead8,R_lead7)= -1;
A1(R_lead9,R_lead9)=1; A0(R_lead9,R_lead8)= -1;
A1(R_lead10,R_lead10)=1; A0(R_lead10,R_lead9)= -1;
A1(R_lead11,R_lead11)=1; A0(R_lead11,R_lead10)= -1;
A1(R_lead12,R_lead12)=1; A0(R_lead12,R_lead11)= -1;


%% Eq. 45-49 Monetary Policy Signals

%% Eq. $$ \zeta^{9}_{t}=\zeta^{10}_{t-1} + \epsilon^{9}_{t} $$ 
A1(sigR_9,sigR_9)=1;PSI(sigR_9,Rs_lag9)=1; 

%% Eq. $$ \zeta^{8}_{t}=\zeta^{9}_{t-1} + \epsilon^{8}_{t} $$ 
A1(sigR_8,sigR_8)=1;A2(sigR_8,sigR_9)=1;PSI(sigR_8,Rs_lag8)=1; 

%% Eq. $$ \zeta^{7}_{t}=\zeta^{8}_{t-1} + \epsilon^{7}_{t} $$ 
A1(sigR_7,sigR_7)=1;A2(sigR_7,sigR_8)=1;PSI(sigR_7,Rs_lag7)=1; 

%% Eq. $$ \zeta^{6}_{t}=\zeta^{7}_{t-1} + \epsilon^{6}_{t} $$ 
A1(sigR_6,sigR_6)=1;A2(sigR_6,sigR_7)=1;PSI(sigR_6,Rs_lag6)=1; 

%% Eq. $$ \zeta^{5}_{t}=\zeta^{6}_{t-1} + \epsilon^{5}_{t} $$ 
A1(sigR_5,sigR_5)=1;A2(sigR_5,sigR_6)=1;PSI(sigR_5,Rs_lag5)=1; 

%% Eq. $$ \zeta^{4}_{t}=\zeta^{5}_{t-1} + \epsilon^{4}_{t} $$ 
A1(sigR_4,sigR_4)=1;A2(sigR_4,sigR_5)=1;PSI(sigR_4,Rs_lag4)=1; 

%% Eq. $$ \zeta^{3}_{t} = \zeta^{4}_{t-1} + \epsilon^{3}_{t} $$ 
A1(sigR_3,sigR_3)=1; A2(sigR_3,sigR_4)=1; PSI(sigR_3,Rs_lag3)=1; 

%% Eq. $$ \zeta^{2}_{t} = \zeta^{3}_{t-1} + \epsilon^{2}_{t} $$ 
A1(sigR_2,sigR_2)=1; A2(sigR_2,sigR_3)=1; PSI(sigR_2,Rs_lag2)=1; 

%% Eq. $$ \zeta^{1}_{t} = \zeta^{2}_{t-1} + \epsilon^{1}_{t} $$ 
A1(sigR_1,sigR_1)=1; A2(sigR_1,sigR_2)=1; PSI(sigR_1,Rs_lag1)=1; 

%% Eq. $$ \zeta^{0}_{t} = \zeta^{1}_{t-1} $$, which implies that the history of
% policy signals regarding the current quarter is given by 
%
% $$ \zeta^{0}_{t} = \sum_{j=1}^{12} \epsilon^{j}_{t-j}  $$
%
% *Note:* The contemporaneous policy innovation appears 
% directly in the definition of the Taylor rule, in case it is AR(1). An
% alternative would be to make the current composite signal +
% contemporaneous innovation persistent. 
A1(sigR_0,sigR_0)=1;A2(sigR_0,sigR_1)=1;  

%% Inflation Expectations Equation
pe_vec = pe_vec(:);
A1(pe_vec,pe_vec)=eye(length(pe_vec)); 
A0(pe_vec,[p;pe_vec(1:end-1)])= -eye(length(pe_vec)); 

A1(expectedPi40,expectedPi40)= -1; 
A1(expectedPi40,pe_vec      )= (expInfLoading/length(pe_vec))*ones(1,length(pe_vec));  

%% Sum of Real Rates Equation
A1(realRate12q,realRate12q)= -1; 
A1(realRate12q,R          )=(1/12)*1; 
A1(realRate12q,Rlead_vec(1:11) )=(1/12)*ones(11,1); 
A1(realRate12q,pe_vec(1:12) )= (1/12)*-ones(12,1);


%% *PART 7: Solve Model using AMASOLVE*
[GG,~,impact,~,~,~,~,eu]=amasolve(A0,A1,A2,zeros(NY,1),PSI) ;
% EU =[1;1] if unique and stable RE solution 
if ~isequal(eu,[1;1]); 
    flag.euok=0; 
    ZZ=[]; C=[]; SDX=[];
    return 
end 

flag.euok=1; 
%% *PART 8: Addition of Measurement Equations for Prices*

%% 8.1 Extend model solution by 2*NSERP rows 
% Number of rows to add to the solution 
nadd=nserp*2; 
% Part 1: Position of prices 
posp =(NY+1):(NY+nserp); 
% Part 2: Position of ID errors 
posid=(posp(end)+1):(NY+2*nserp); 
GG=blkdiag(GG,zeros(nadd)); 
% Matrix of impact coefficients 
RR=zeros(NY+nadd,NX+nserp); 
RR(1:NY,1:NX)=impact; 
impact=[]; 
%% 8.2 Idiosyncratic errors 
% Which evolve as $$ \nu^{i}_{t} = \rho  \nu^{i}_{t-1} + \varsigma^{i}_{t} $$
% as independent AR(1) processes 
GG(posid,posid   )=diag(rhoid); 
RR(posid,NX+1:end)=eye(nserp); 

%% 8.3 General Structure of Observable Price Equation 
% 
% For price $$ P^{i,C}_{t} = \beta^{i} P^{C}_{t} + \nu^{i}_{t} $$ written in terms
% of lagged variables as 
% $$ P^{i,C}_{t} = \beta^{i} G_{P^{C}} S_{t-1} + \beta^{i} R_{P^{C}}
% \eta_{t} + \rho \nu^{i}_{t-1} + \varsigma^{i}_{t} $$ ,
% where $$ \beta^{i} $$ is the series specific factor loading, and, $$ G_{P^{C}}  R_{P^{C}} $$ correspond to the rows of the solution
% matrix for the model consumption price 

%% 8.3.a First row is GDP 
GG(posp(1),1:NY)=GG(gdpdef,1:NY);
RR(posp(1),1:NX)=RR(gdpdef,1:NX); 

%% 8.3.b Remaining rows are multiple indicators for consumption prices 
% Recall Structural part of observation equation $$ \beta^{i} G_{P^{C}} S_{t-1} +
% \beta^{i} R_{P^{C}} \eta_{t}  $$ 
GG(posp(2:end),:)       =repmat( GG(p,:)    ,[nserp-1 1] ); 
RR(posp(2:end),1:NX)    =repmat( RR(p,1:NX) ,[nserp-1 1] );
% Multiply by Factor Loadings  
bet_c=bet_c(:); 
GG(posp(2:end),1:NY)=(repmat(bet_c,[1 NY])).*GG(posp(2:end),1:NY); 
RR(posp(2:end),1:NX)=(repmat(bet_c,[1 NX])).*RR(posp(2:end),1:NX);
%% 8.3.c Add idiosyncratic errors 
GG(posp,posid)   =diag(rhoid); 
RR(posp,NX+1:end)=eye(nserp); 

% Update model dimensions after adding prices 
%impact=RR; 
NY=size(GG,1); 
NX=size(RR,2); 

%% *Part 9: Addition of Measurement Equations for Spread*
% 
%% 9.1 Extend model solution by 1 row 
GG    =blkdiag(GG,zeros(2));
RR   =[RR zeros(NY,1);zeros(2,NX+1)];

%% 9.2 Fill in last row for the spread 
% 
% $$ sp_{t} = \kappa \upsilon_{t} + \varsigma^{sp}_{t} $$ 
posspread=NY+1;
GG(posspread,:)   =kappa*GG(upsilon,:); 
RR(posspread,:)   =kappa*RR(upsilon,:); 
RR(posspread,NX+1)=1; 

%% 9.3 Fill in spread measurement error
% 
RR(posspread+1,NX+1)=1; 

% Update model dimensions after adding spread
NY=NY+2;
NX=NX+1;


%% *Part 10: Declare Real Observable {Y,C,I,W,SumR}, Accounting for Measurement Errors* 
nadd2=4; 
GG       =blkdiag(GG,zeros(nadd2)); 
RR       =[RR;zeros(nadd2,size(RR,2))]; 
posGDPreal  =NY+1; 
posCreal    =NY+2; 
posIreal    =NY+3; 
posWreal    =NY+4; 
%posSumR    =NY+5;

%% Eq. Real GDP: Nominal Deflated by GDP deflator 
GG(posGDPreal,:)     =GG(ynomobs,: )-GG( posp(1),:); 
RR(posGDPreal,:)     =RR(ynomobs,: )-RR( posp(1),:); 

%% Eq. Real Consumption: Nominal Deflated by C deflator 
GG(posCreal,:) =GG(cnomobs,: )-GG( posp(3),:); 
RR(posCreal,:) =RR(cnomobs,: )-RR( posp(3),:); 

%% Eq. Real Investment: Nominal Deflated by I deflator 
GG(posIreal,:) =GG(inomobs,: )-GG(dpinv,:); 
RR(posIreal,:) =RR(inomobs,:) -RR(dpinv,:); 

%% Eq. Real Wages: Nominal Deflated by C deflator 
GG(posWreal,:) =GG(wnomobs,: )-GG(posp(3),:); 
RR(posWreal,:) =RR(wnomobs,: )-RR(posp(3),:); 

% %% Sum of Expected Real Rates over 3 Year Horizon
% GG(posSumR,:) = GG(R_lead1,:)+GG(R_lead2,:)+GG(R_lead3,:)+GG(R_lead4,:)+GG(R_lead5,:)+GG(R_lead6,:)+...
%                 GG(R_lead7,:)+GG(R_lead8,:)+GG(R_lead9,:)+GG(R_lead10,:)+GG(R_lead11,:)+GG(R_lead12,:)-...
%                 GG(pe_vec(2),:)-GG(pe_vec(3),:)-GG(pe_vec(4),:)-GG(pe_vec(5),:)-GG(pe_vec(6),:)-GG(pe_vec(7),:)-GG(pe_vec(8),:)-...
%                 GG(pe_vec(9),:)-GG(pe_vec(10),:)-GG(pe_vec(11),:)-GG(pe_vec(12),:)-GG(pe_vec(13),:);  
NY=size(GG,1); 

%% *PART 11: Declare Cholesky of Variance Covariance Matrix* 
% Recall $$ V = SDX'SDX $$ and the order of the shocks is 
% 1. Structural 2. Policy Signals 3. Idiosyncratic Prices 4. Measurement
% Error Spread 
SDX=diag([sdstruct(:);stdsignals(:);sdid(:);sdspreadme]);
if size(SDX,1)~=NX;error('Dimensions of errors and SDX do not match');end

%% *PART 11: Declare ZZ, Matrix of Pointers to Observables*
% Ordering 
% 1) Nominal GDP Growth. 
% 2) Nominal C   Growth.
% 3) Nominal I   Growth.
% 4) Hours. 
% 5) Nominal W   Growth. 
% 6) Feds Fund Rate. 
% 7) Investment Deflator. 
% 8) GDP Deflator. 
% 9 though NSERP-1) Indicators of C Prices. 
% Last is Spread
dimZ=9+nserp+NExFFRObserved; 
ZZ=zeros(dimZ,NY); 
ZZ(1,ynomobs)=1; 
ZZ(2,cnomobs)=1; 
ZZ(3,inomobs)=1; 
ZZ(4,L)=1; 
ZZ(5,wnomobs)=1; 
ZZ(6,R)=1;
ZZ(7,dpinv)=1;
endZPrices=8+nserp-1; 
ZZ(8:endZPrices,posp)=eye(nserp);
ZZ(endZPrices+1,posspread)=1; 
ZZ(endZPrices+2,expectedPi40)=1; 
if endZPrices+2+NExFFRObserved ~=dimZ; error('Wrong Z matrix'); end 

%% Number of Expected FFR that are observable 
if NExFFRObserved > 0
    %Recall Rlead_vec stored the position of ExFFR1 through ExFFR12
    endZnonExFFR=endZPrices+2;
    posZExFFR   =(endZnonExFFR+1:1:endZnonExFFR+NExFFRObserved);
    ZZ( posZExFFR , Rlead_vec(1:NExFFRObserved) )=eye( NExFFRObserved );
end

%% *PART 12: Declare C, Matrix of Constants*
C=zeros(NY,1); 
C(ynomobs)=gammastar100+pss100;
C(cnomobs)=gammastar100+pss100;
C(inomobs)=gammastar100+pss100;
C(L)=Lss;
C(wnomobs)=gammastar100+pss100;
C(R)    =pss100+rss100;
C(dpinv)=(pss100-gammamiu100)+cons_idef; 
C(posp )= pss100*ones(nserp,1)+cprices(:); 
C(expectedPi40)=pss100+expInfCons ;
if NExFFRObserved > 0
    C(Rlead_vec(1:NExFFRObserved) )=C(R)*ones(NExFFRObserved,1);
end
if flag_repstanames==0; 
    return 
end 
%________________________________________________________
%% 16.1 STANAMES: Cell with Names of ALL states 
stanames1={'yhat','khat','hours','Return K','what','x','dp','Marginal Cost','lambdahat','chat',...
    'nomr','Utilization ','phihat','ihat','kbar','mp','z','g','miu','lambdap',...
    'psi','b','upsilon ','pitarg','wmarkwn','GDPm nominal','Cm nominal','Im nominal','Wm nominal','GDPm real',...
    'GDPm level','I deflator','GDPm deflator','Cm real','Im real','Wm real','Output gap','Real rate','GDPm lag1','GDPm lag2',...
    'GDP annualized','Inflation lag1','Inflation lag2','Inflation annualized','Aggregate trend','ExFFR1','ExFFR2','ExFFR3','ExFFR4',...
    'ExFFR5','ExFFR6','ExFFR7','ExFFR8','ExFFR9','ExFFR10','ExFFR11','ExFFR12','Signal0','Signal1','Signal2',...
    'Signal3','Signal4','Signal5','Signal6','Signal7','Signal8','Signal9','Signal10','Signal11','Signal12',...
    'ExpInf1','ExpInf2','ExpInf3','ExpInf4','ExpInf5','ExpInf6','ExpInf7','ExpInf8','ExpInf9','ExpInf10',...
    'ExpInf11','ExpInf12','ExpInf13','ExpInf14','ExpInf15','ExpInf16','ExpInf17','ExpInf18','ExpInf19','ExpInf20',...
    'ExpInf21','ExpInf22','ExpInf23','ExpInf24','ExpInf25','ExpInf26','ExpInf27','ExpInf28','ExpInf29','ExpInf30',...
    'ExpInf31','ExpInf32','ExpInf33','ExpInf34','ExpInf35','ExpInf36','ExpInf37','ExpInf38','ExpInf39','ExpInf40','AvgExpInf10y','RealRate12q',...
    'ystar','kstar','hstar','rkstar','wstar','xstar','rstar','mgcstar','lambdastar','cstar',...
    'ustar','phistar','istar','kbarstar','GDP star','GDP deflator'}; 
if nserp==3 
    stanames2={'C deflator','PCE Core','GDP def Id','C def Id','PCE Core Id'}; 
elseif nserp==4 
    stanames2={'C deflator','PCE Core','CPI Core','GDP def Id','C def Id','PCE Core Id','CPI Core Id'}; 
else 
    disp('Nprices > 4 please add names to STANAMES') 
end 
stanames3={'Spread','Spread Me','RealY','RealC','RealI','RealW'}; 
stanames=[stanames1(:);stanames2(:);stanames3(:)]; 

%% 16.2 STSHOCKS: Cell with Names of Shocks in STANAMES 
% *Note:* Technically not true for the policy signals, for which looking
% only at the cummulative signals. Hence isolate policy innovations rather
% than policy states.
stshocks={'mp','z','g','miu','lambdap','psi','b','upsilon ','pitarg','wmarkwn',...
    'Signal1','Signal2','Signal3','Signal4','Signal5','Signal6','Signal7','Signal8','Signal9','Signal10','Signal11','Signal12',...
    'GDP def Id','C def Id','PCE Core Id'};
if nserp==4
    stshocks=[stshocks,'CPI Id'];
end
stshocks=[stshocks,'Spread Me'];
stshocks=stshocks(:);

%% 16.3 SHONAMES: Cell with Names of Shock Innovations 
shonames={'MP','Neutral','G+NX','ISTS','Price Markup','Labor Disutility',...
    'Discount','MEI','Inflation Drift','Wage Markup',...
    'MP Signal 1','MP Signal 2','MP Signal 3','MP Signal 4','MP Signal 5','MP Signal 6','MP Signal 7',...
    'MP Signal 8','MP Signal 9','MP Signal 10','MP Signal 11','MP Signal 12'}; 
shonames=[shonames(:);stshocks(23:end)]; 


%obsnames={}; for ii=1:size(ZZ,1); temp=stanames( find(ZZ(ii,:)~=0 ) ); obsnames=[obsnames;temp(:)]; end;

%% 16.4 Quick check dimensions 
if length(stanames)~=size(GG,1)
    error('Mistmatch STANAMES and NY');
end
if length(stshocks)~=size(SDX,1)
    error('Mistmatch STSHOCKS and NX');
end
if length(shonames)~=size(SDX,1)
    error('Mistmatch SHONAMES and NX');
end
